A&A manuscript no. 

(will be inserted by hand later) 



Your thesaurus codes are: 

02.01.1, 11.01.2, 11.14.1, 11.19.1, 13.07.3, 13.25.2 



ASTRONOMY 
AND 
ASTROPHYSICS 
19.9.1994 



Self- Consistent Particle Acceleration in Active Galactic 
Nuclei 

A. Mastichiadis and J. G. Kirk 

Max-Planck-Institut fiir Kernphysik, Postfach 10 39 80, 69029 Heidelberg, Germany 
Received 15.7.94, accepted 7.9.94 



Abstract. Adopting the hypothesis that the nonthermal emis- 
sion of Active Galactic Nuclei (AGN) is primarily due to the 
acceleration of protons, we construct a simple model in which 
the interplay of acceleration and losses can be studied together 
with the formation of the emitted spectrum. The accelera- 
^— H tion process is assumed to be of the first order Fermi type, 
^ and the proton distribution as well as the injected electrons 
^^and photons in the central region of the AGN are described 
by spatially averaged kinetic equations. The various relevant 
processes which dominate the three species are incorporated 
into the equations. The technique used to solve these is pre- 
' sented and several tests of the numerical implementation are 
presented. We also present results of a sample time-dependent 
AGN model in which photons appear suddenly as a result of a 
— feedback instability and the system evolves to a steady state, 
^ 1 ^ in which the acceleration process is saturated self-consistently 
00 by the photons it produces. This example combines an X-ray 
^ power law index of about —1.7, together with a break at an 
^ energy between 50 and 500 keV. 
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1. Introduction 

Active Galactic Nuclei are thought to be powered by accretion 
of matter onto a central massive black hole (Rees 1984). How- 
ever, the way in which gravitational energy is converted into 
electromagnetic radiation remains largely unclear. Studies of 
the X-ray and 7-ray properties of AGNs can be particularly il- 
luminating as variability (especially in the well-measured X-ray 
regime) implies that this radiation is generated within a few 
Schwarzschild radii of the black hole (Done &; Fabian 1989). 

Over the last decade, considerable effort has been expended 
on understanding the X-ray spectrum of AGNs. The high val- 
ues which X-ray observations imply for the photon compact- 
ness ij = i^(TT/4iriZmeC^ (where Lj is the luminosity, R 
the radius of the source and ctt the Thomson cross-section) 
mean that any 7-ray of energy > 1 MeV, if indeed present, 
would be readily absorbed by the softer photons. The result- 
ing electron-positron pairs would quickly lose their energy, pro- 
ducing more energetic photons and an intense electromagnetic 
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cascades would ensue (for a review see Svensson 1989 and ref- 
erences therein). Such a cascade could be expected to alter 
drastically the primary photon spectrum, and it was hoped 
that this property might explain the canonical X-ray power 
law spectral index of —1.7, as observed by both the HEAO-1 
(Rothschild et al. 1983) and the GINGA (Turner fc Pounds 
1989) X-ray observatories. However, these models do not in- 
clude a discussion of particle acceleration, but instead describe 
the required injection of energetic electrons and photons by 
a number of free parameters. It turns out that only a rather 
limited region in this parameter space gives acceptable results 
(Lightman &; Zdziarski 1987) - a conclusion which is, however, 
relaxed by the introduction of models including reflection (see 
Pounds et al. 1989 and Zdziarski et al. 1990). Another short- 
coming of pair cascade models is that although variability is 
one of the basic characteristics of AGNs most of them assume 
steady state conditions. Without a model of acceleration, it is 
possible to investigate only rather artificial variations of the 
injection paratmeters (Done &; Fabian 1989, Coppi 1992). 

Approaching the problem from a different angle, a class 
of models proposing the presence of a strong stationary shock 
in the central part of an AGN has been developed (Protheroe 
&; Kazanas 1983, Kazanas &; Ellison 1986). Such a shock can, 
in principle, accelerate protons to very high energies, thereby 
converting part of the gravitational energy released by the 
plasma accreting into non-thermal particles. The presence of 
high energy protons in AGNs has many interesting conse- 
quences, which have been investigated in a number of pa- 
pers: relativistic protons will eventually lose part of their en- 
ergy in inelastic collisions with ambient protons (Protheroe &; 
Kazanas 1983, Kazanas &; Ellison 1986) or ambient photons 
(Sikora et al. 1987). In each case neutrons and neutrinos, as 
well as high energy electrons/positrons and 7-rays will be pro- 
duced. Unlike protons, which are magnetically confined inside 
the acceleration region, the neutrons escape, carrying off a frac- 
tion of the initial luminosity (Kirk &; Mastichiadis 1989, Sikora 
et al. 1989, Atoyan 1992a). Various observable or potentially 
observable effects of these neutrons have been predicted, such 
as the production of Boron in spallation reactions (Kirk &; 
Mastichiadis 1989), the emission of Very High Energy 7-rays 
(Mastichiadis &; Protheroe 1990, Atoyan 1992b), the acceler- 
ation of winds in broad absorption line QSOs (Begelman et 
al. 1991) and the production of cosmic rays around the 'knee' 
(Protheroe & Szabo 1992). Also, the flat radio spectra of ra- 
dio loud AGNs have been attributed indirectly to escaping 



2 



neutrons (i.e., after their decay into protons and electrons - 
Giovanoni &; Kazanas 1990). Finally, the flux of escaping neu- 
trons has been used to calculate the contribution of AGNs to 
the diffuse 7-ray background (Johnson et al. 1994). Neutrinos, 
on the other hand, may be produced at a rate which could 
be detected with an experiment such as DUMAND (Stecker 
et al. 1991, Biermann 1992, Sikora &; Begelman 1992, Szabo &; 
Protheroe 1992, 1994). In these models the relativistic protons 
also inject a population of electrons and positrons, which cool 
by synchrotron and/or inverse Compton radiation and initiate 
intense pair cascades. However, the resulting nonthermal radi- 
ation (which is, at least partly, responsible for the saturation of 
the acceleration) was specified a priori and not calculated self- 
consistently. In this sense, the hadronic models include a model 
of particle acceleration and injection but lack a self-consistent 
treatment of the electromagnetic cascades. 

A first attempt to create a synthesis of these models was 
made by Stern et al. (1991), Stern &; Svensson (1991) and Stern 
et al. (1992). These authors followed the electromagnetic cas- 
cades resulting from the injection of electrons by relativistic 
protons and included the feedback of the photons on the rel- 
ativistic protons. They found that this system showed limit 
cycles much like a predator-prey system. However, the Monte 
Carlo approach which they use complicates the interpreta- 
tion of the results in terms of a specific feedback mechanism. 
Such a feedback effect was found analytically by Kirk &; Mas- 
tichiadis (1992 - henceforth "KM92") using the kinetic equa- 
tion approach. They showed that relativistic protons become 
unstable to a combination of proton-photon pair production 
and electron synchrotron radiation once their number density 
and energy exceeds a certain critical value. This feedback leads 
eventually to rapid energy losses for the relativistic protons. 
However, the analysis of KM92 employs a stationary proton 
distribution and uses a linearised set of equations which ex- 
clude electromagnetic cascades, making it impossible to follow 
the evolution of the system into the saturation phase. 

Motivated by this shortcoming we present here a model ca- 
pable of describing time-dependent effects in the production of 
the nonthermal spectra of AGNs, albeit under a set of highly 
simplifying assumptions. Our method is to describe the three 
basic components of the central region of an AGN - protons, 
electrons and photons - by a system of three spatially aver- 
aged kinetic equations. The aim is to incorporate in an ap- 
proximate manner all the important processes acting on these 
constituents and to use numerical methods to integrate the 
system forwards in time in the manner described by Fabian 
et al. (1986) and Coppi (1992) for photons and electrons. A 
time-dependent method is also essential if one is to account 
for the highly nonlinear coupling of the acceleration process 
with losses caused by the associated photons. We propose to 
use this technique to gain a better understanding of the ori- 
gin and properties of variability in AGNs as well as the way 
in which their photon spectra are formed. Our method should 
also provide better information about quantities such as the 
expected luminosity in high energy neutrinos. In the present 
paper, however, we restrict ourselves to a detailed description 
of the method together with a discussion of its strengths and 
weaknesses. We do not attempt a systematic investigation of 
the parameter regime appropriate to AGNs, but present a sam- 
ple set of results obtained using the full code. Preliminary 
results of this work have been presented by Mastichiadis &; 
Kirk (1992). 



The contents of the paper are organised as follows: in 
Sect. 2 we describe the way in which we model the acceleration 
process. Only protons are assumed to undergo acceleration - it 
being tacitly assumed that the relativistic electrons are dom- 
inated by rapid loss processes. We choose a model in which a 
first-order partial differential equation is used to describe the 
first-order Fermi process and a 'loss' or 'escape' probability is 
introduced to account for the possibility that protons might 
leave the emission region, for example by accretion into the 
black hole. In such a model, acceleration occurs homogeneously 
throughout the emission region, as might be expected, for ex- 
ample, if a converging accretion flow provides the acceleration 
(Schneider &; Bogdan 1989). 

Apart from acceleration, the microscopic processes of im- 
portance in the system of kinetic equations are relatively well 
understood. For the protons these include proton-proton and 
proton-photon collisions (Mannheim &; Biermann 1989, Begel- 
man et al. 1990) whereas the electrons and photons (in ad- 
dition to the source terms provided by the proton related 
processes) experience synchrotron radiation, Compton scatter- 
ing, photon-photon pair production, electron-positron annihi- 
lation and Compton downscattering on cooled electrons (Coppi 
&; Blandford 1990). These processes and the approximations 
we employ to describe them are discussed in Sect. 3. 

The method used to convert the integro-differential kinetic 
equations into a system of ordinary differential equations suit- 
able for integration by standard numerical methods is pre- 
sented in Sect. 4, and this is followed by a series of tests which 
check the behaviour of our approximation schemes in circum- 
stances in which either analytic solutions are available e.g., 
when only synchrotron cooling or inverse Compton cooling 
are present, or in which there exist calculations in the liter- 
ature with which to compare, e.g., the spectra of stationary 
electromagnetic cascades in the Thomson regime (Lightman &; 
Zdziarski 1987). 

As an example of the application of the full code. Sect. 5 
presents time-dependent results for one particular set of pa- 
rameters which are appropriate for an AGN. One of the fea- 
tures of this run is the development of the pair production- 
synchrotron instability (KM92) once the marginal stability 
threshold is crossed. The X-ray spectrum of power-law index 
a = —5/3 which is predicted in the linear phase of the insta- 
bility is seen to persevere well into the nonlinear phase and 
a stationary state is achieved in which the photons produced 
by accelerated protons are responsible for saturating the accel- 
eration process. To conclude, we briefly summarise the main 
advantages and limitations inherent in our approach in Sect. 6. 

2. The Particle Acceleration Model 

According to the theory of particle acceleration by the first- 
order Fermi mechanism (see, for example. Kirk et al. 1994) 
stochastic encounters of particles with so-called 'scattering cen- 
tres' occur within an 'acceleration region', such that each in- 
teraction results in a small increase in the particle's energy. 
At the same time, a particle has a finite probability of es- 
caping from the acceleration region, and the combination of 
these two effects leads to a distribution of accelerated particles 
which, under certain circumstances, is of the power-law type. 
We shall adopt this theory here, and make the simplest possi- 
ble assumptions concerning the rate at which energy is gained 
from the scattering centres as well as the rate at which parti- 
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cles escape the region containing these centres. The aim is to 
describe the particle distribution using a kinetic equation in 
which not just acceleration, but also losses suffered as a result 
of interactions with ambient photons can be included. However, 
one important question about the nature of an AGN, which we 
cannot answer a priori is that of how large the acceleration re- 
gion is compared to the size of the 'source region', i.e., that 
region which contains the relativistic particles responsible for 
the observed nonthermal emission. Two possibilities represent 
opposite extremes: 

1. acceleration occurs throughout the whole of the source re- 
gion, and 

2. acceleration takes place in only very small parts of the 
source. 

An example of case 1 is the acceleration of particles by a 
smooth, unshocked accretion flow (Payne &; Blandford 1981, 
Cowsik fc Lee 1982, Webb fc Bogdan 1987, Schneider fc Bog- 
dan 1989). Particles may escape from the acceleration region 
in this case (perhaps by being accreted into the black hole) 
but if they do, they cease to emit observable radiation and, 
by definition, have left the source. An example of case 2 is the 
acceleration of particles at a shock front in the accretion flow 
(Protheroe &; Kazanas 1983, Sikora et al. 1987). Provided the 
mean free path of the particle is small compared to the source, 
only those particles in the immediate vicinity of the shock un- 
dergo acceleration. Once particles have been swept out of this 
region, they are highly unlikely to return to it, but may still be 
energetic enough to produce observable radiation whilst cool- 
ing on the ambient photons, and so cannot be considered to 
have left the source. In this paper we deal exclusively with 
case 1, in which scattering centres are distributed throughout 
the source region. 

Denoting by jip(p)dp the differential number density of pro- 
tons of momentum p in the interval dp, the kinetic equation for 
protons within the source region can be written (e.g., Schlick- 
eiser 1984, Kirk et al. 1994) 



dt dp 



- — np{p,t) 

tacc 



iesc 



Qinj*(p -Pinj) + C^{np,p,t) 



(1) 



where Qinj is the number of particles injected into the acceler- 
ation process per second per unit volume with momentum pinj 
and (which can be a differential and/or integral operator 
acting on rip) denotes the losses suffered by energetic protons. 
The second term in Eq. (1) provides a continuous energy in- 
put by the first-order Fermi process into those protons which 
remain in the acceleration region, whereas the third allows for 
escape, the average residence time being t^sc- 

If we ignore for a moment the losses, and take the sim- 
ple case of constant Qinj, tacc, and tesc, the solution to Eq. (1) 
which satisfies the boundary condition jip(p, 0) = is (Ax- 
ford 1981) 



np{p,t) 



^acc^inj / P 



-1-t.oo/tc 



[if(p-p;nj)-if(p-p;nje'/'>")] 



(2) 



where H[x) is the Heaviside function, equal to zero for a; < 
and unity for a; > 0. This solution is just a power-law extending 



from Pinj up to a cut-off at pmax(t) = Pinje''''""" which increases 
with time. Below the cut-off i.e., for p < pmax(t) the solution 
is independent of time. The power-law index is determined by 
the relative strengths of the acceleration term and escape term. 

In the following we will usually assume tacc = tesc, in which 
case hp oc p~^ , such as is expected, for example, of particles 
accelerated at a strong shock front in a gas of adiabatic in- 
dex 5/3. Interesting conclusions can be drawn in this special 
case about the energy given to and extracted from the pro- 
tons by building the moment of Eq. (1) with the kinetic en- 
ergy of a proton: (7 — l)mpC , where 7 is the Lorentz factor 
(7 = -^1 +p^/(mpc)^). Denoting the total energy (minus the 
rest mass) contained in accelerated protons by E, we find after 
integrating by parts: 



dE 
dt 



FQi„jmpC^(7i„j - 1) - iioss 



-j-VnipC' 
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dp 



(7 - l)hp{p,t) 



(3) 



where 7inj = ^1 +Phij/("*p'-)^' ^loss is the rate at which en- 
ergy is lost by protons (to the processes of pair production 
and pion production) and V is the source volume. The loss 
processes discussed in the following section operate effectively 
only on relativistic protons, so we can expect that if a steady 
state is set up in which E = constant, the spectrum will be 
given by the loss-free solution Eq. (2) all the way from the in- 
jection momentum up into the relativistic regime. However, the 
integral in the third term on the right-hand side of Eq. (3) is 
dominated by the nonrelativistic and transrelativistic regimes, 
provided the particle density does not diverge at large p. Con- 
sequently, we make only a small error by using the loss-free 
distribution in this integral. Provided pinj <C nipC, we find in 
the steady state: 



iloss = l^CpinjQinj. 



(4) 



According to this equation, the rate at which protons put en- 
ergy into pair production and pion production during the ac- 
celeration process (which is in this model the entire nonther- 
mal luminosity of the AGN) is determined in the steady state 
solely by the rate at which they are injected into the accel- 
eration process at low momentum. In particular, it is inde- 
pendent of quantities connected with the actual loss process 
itself, such as the background photon or matter density and 
the maximum Lorentz factor to which particles can be acceler- 
ated, even though these may depend nonlinearly on the density 
itself. In connection with Eq. (4), we note that the rate iinj at 
which energy is injected at momentum pinj is small compared 
to the nonthermal luminosity: iinj/iioss = Pinj/2mpC <C 1, so 
the nonthermal emission stems not from the unknown injec- 
tion process, but from the first-order Fermi mechanism we are 
modelling. 

In this study of AGNs, we will be concerned only with 
relativistic protons, since it is in this regime that the accelera- 
tion and loss processes can compete with each other. It is then 
more convenient to write the kinetic equation (1) in terms of 
the Lorentz factor 7 of the protons, using a normalisation in 
which time is measured in units of the light crossing time of 
the source (of size R) and the particle density refers to the 
number contained in a volume element of size (TtR (where ctt 
is the Thomson cross section). Accordingly, we define the di- 
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mensionless density np(7,t) by 

np{~/,t)d~/ = (TT-Rjip(p, t)dp 

and find, in the relativistic regime: 



(5) 



dt 



iesc 



+ ^ [jt^^i^'')] + '-"^7^ = 'C-K.T,*). (6) 



where t is now dimensionless, tacc = ctg.cc/R and tcsc = ctcsc/R- 
The nonrelativistic part of the acceleration process can be 
avoided by using the loss-free solution as a boundary condi- 
tion 



/ .\ / l + tacc/tcsc 

"p(7o,*) = "o/To I 
with 

no = <TT-R4ccQinj(pinj/mpC)'""'''°" , 



(7) 



(8) 



where 70 is the lowest value of the Lorentz factor to be con- 
sidered. For tacc = tcsc, the AGN luminosity in a steady state 
can be expressed in terms of the compactness 



4ir iZm.c^ 



Setting the source volume V = jT^R^ in Eq. (4) leads to 
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no f nip \ 
3tacc \mcJ 



(9) 



(10) 



Using the same normalisation, we complement Eq. (6) by 
writing for the relativistic electrons and positrons: 



dnc{~f,t) 
dt 



Q^{nc,~f,t) + C{nc,~f,t) 



Here C denotes the various electron loss terms while Q' are 
the injection terms. Note that there is no acceleration or escape 
term included in the electron equation, since we assume the loss 
terms to be much larger. 

In our numerical treatment, we impose a lower limit 7min 
on the Lorentz factor of the relativistic electron population and 
assume particles which cool through this boundary join a pop- 
ulation of cold electrons whose number density Nl°°^(t) (i.e., 
the number in a volume ctt-R) is determined by the equation 



J Arcool/j\ 

^-^'e \^/ ^e,cool/ .\ 1 /»e,cool/ .\ 

= Q {nc,t) + C (nc,t) 



(12) 



Contributions to the source term Q°''^°°'(ne, t) arise from syn- 
chrotron and Compton cooling, while the contributions to the 
sink term £°''^°°'(ne, t) come mainly from electron-positron an- 
nihilation. 

Finally, the spatially averaged photon equation reads: 

dn^(x,t) n^(x,t) . „„, . , . 

"2, + I = Q^n^,x,t) + C^n^,x,t) , (13) 

where x is the dimensionless photon frequency: x = hu/[mcC^). 
Photons leave the source on the timescale t-ycsc (measured in 
units of the light crossing time tcross = R/c) and, rather than 
being advected into the black hole, are assumed to propagate 
freely after escape. The terms C and Q'' denote the sinks 



and sources of photons. Because of the normalisation used, the 
quantity J dx xn^ is simply related to the photon compactness 
(Guilbert et al. 1983) given by 
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For a spherical source we have 
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3. The Physical Processes 

We proceed now to discuss the various contributions to the ki- 
netic equations (6), (11), (12) and (13) of the physical processes 
by which the components of the system interact with each other 
and with the magnetic field. 

3.1. Proton-proton interactions 

Inelastic proton-proton collisions act as an energy loss mech- 
anism for relativistic protons. They also inject 7-rays, rela- 
tivistic electrons, and neutrinos resulting from the decay of 
the produced neutral and charged pions. While one can cal- 
culate in detail the spectra of the products once the rela- 
tivistic proton distribution is given (e.g., Dermer 1986, Mas- 
tichiadis &; Protheroe 1990), for the present calculation it suf- 
fices to use the 5— functional approximation to the differential 
cross-section for production of energetic pions da^i^E-^) I dE-^ = 
(Tpp(TT5(-Bir — 0.15 7mpC ) where Upp ~ 0.06 is the proton- 
proton cross section (e.g., Atoyan 1992a). 
Thus, to treat the proton losses we write 



(11) ^pp(7, *) = <^pp"p"* ["pCt', *) - "p(7, *) 



(16) 



where nj,*'* is the density of target protons and we have as- 
sumed that each proton-proton collision removes protons from 
the energy bin 7 and 7 + ^7 but injects them there from higher 
energy bins of energy 7' = 7/(1 — fcpp). Here fcpp is the proton 
inelasticity and it is taken to be fcpp = .45 

Since the mean energy per gamma-photon produced at the 
decay of the ir^-meson is O.^Et, = 0.0757mpC^, the photon 
production spectrum, in this approximation, can be taken as 



^{x,t) 



0.075 



"p"*"p(7i,*) 



(17) 



where 71 = a;me/(0.075mp). The above approximation essen- 
tially injects photons with a slope equal to the relativistic pro- 
ton spectrum at energies greater than ~ 70 MeV. It does not 
treat correctly injection at energies below this value. However, 
detailed calculations (e.g., Dermer 1986) have shown that the 
photon spectrum flattens considerably at energies below 100 
MeV, so that the approximation introduced above is adequate 
for our purposes. 

We turn next to the injected electrons (or positrons). Since 
the electron (positron) carries off, on average, 26% of the initial 
pion energy, it follows that {Ec) ~ 0.039 E-p (Atoyan 1992a) 
and the corresponding injection spectrum resulting from the 
decay of charged pions is given by 



2;p(7,<) 



0.039 



tare / ,\ 

CTppiip ^np{~f2,t) 



(18) 
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where 72 = 7me/(0.039mp). 

Using similar arguments one can calculate the neutrino 
emissivity resulting from proton-proton collisions. Thus, we 
write 



0.037 



<Tppnp"Snp(73,t) , 



(19) 



where 73 = £;„/(0.037mpC^). 

It is easy to show using the above relations that the energy 
lost per unit time by protons is equal to the amount of energy 
per unit time given to photons, electrons and neutrinos, i.e., 
our approximate treatment conserves energy. 

Finally, it is interesting to note that if we assume the proton 
losses in Eq. (16) are catastrophic, i.e., if we write 



£pp(7,*) = -<^pp"p"* ["p(7,*)] 



(20) 



then instead of Eq. (16), we obtain essentially Eq. (6) and can 
immediately write the analytic solution: 



"p(7,*) = "0 ( ^ 



l + (t.oo/tcso) + «p"*<T° t.oo 



[-ff(7-7inj)--ff(7-7inje'/'"") 



(21) 



Comparing the above solution with Eq. (2) we see that the 
inclusion of proton-proton losses causes the proton distribution 
to become steeper, but does not affect the time evolution of the 
upper cut-off. 

3.2. Proton-photon pair production 

A photon of (dimensionless) energy x can produce an elec- 
tron/positron pair in the Coulomb field of the proton if the 
threshold condition 



7pa; > 2 



(22) 



is satisfied. Assuming that the resulting electron and positron 
have the same Lorentz factor as the incoming proton (which 
is true provided the proton-photon collisions occur predom- 
inantly close to threshold), the fractional energy loss of the 
proton is small. In this case the losses can be considered a 
continuous process and can be written as 



£?^^pee(7, t) = 2^ A (7rpnp(7, t)) 



where 



P oc 



df 



(23) 



(24) 



is the collision rate and <Tpe(j/) is the cross-section of the process 
in units of the Thomson cross-section as a function of the pho- 
ton energy y as seen in the rest frame of the proton. 

On the other hand, the rate at which this process in- 
jects electrons and positrons, (we make no distinction between 
them) is 



3.3. Proton-photon pion production 

To take this complicated process into account we assumed that 
the cross section is given by a 5-function, i.e., d(Tp„(j/)/dj/ = 
''■pir'''T5(j/ — yo), where again y is the photon energy seen from 
the proton rest frame, (Tp„ = .25 and j/o = 10 . We also take 
the inelasticity to be fcp = .2 and the multiplicity to be 1. 
These values are found to be in good agreement (better than 
10%) with the results of Begelman et al. (1990) who calculate 
the proton losses due to photopion production using the full 
cross section. We consider the two basic channels 



(a) 



(b) 



p + 7 ^ p + ir 



p + 7 ^ n + ir 



(26) 



(27) 



to be equally probable. Whereas the outgoing protons of chan- 
nel (a) are effectively trapped, the neutrons of channel (b) 
are assumed to escape the source without further attenuation. 
Thus we treat channel (a) as an energy loss process which 
preserves proton number in contrast to channel (b) which is 
treated as a catastrophic proton loss. The photopion collisions 
of channel (a) move protons of Lorentz factors between 7 and 
7 + d7 to lower ones and add photons to this range which, be- 
fore interaction, had Lorentz factors around 7' = 7/(1 — fcp). 
Thus, using the 5-function approximation for the cross section, 
we find for channel (a): 



■^P7— >p7r,a (7? 



-p7->p 

2 



^np{-y',t)n^{yo/nf',t) - -np{-y,t)n^{yo/nf,t) 
7 7 



and for the catastrophic losses of channel (b): 



P7 



^p,r,b(7i*) = ^-np{nf,t)n^{yo/nf,t) 



(28) 



(29) 



As neutral pions decay essentially instantaneously into 7- 
rays, channel (a) will provide a source term in the photon equa- 
tion Eq. (13). Assuming the two 7-rays have equal energy (see, 
for example, Stecker 1968), this quantity (a;) is related to the 
incoming proton energy 71 by a; = (mp/2me)fcp7i . Thus, we 
can write the photon source term as 



Qpi^p-„{x,t) = np{-yi)n^{xi)-^ (30) 

where xi = kpUipyo /{2xmc). 

Three neutrinos and a positron are created in the decay 
chain of a ir'*' from channel (b). Assuming again that these are 
produced with equal energies, we can write the electron source 
term analogous to the photon term in Eq. (30): 



Qp'y—^p-jT (7? 



1 f 

-"p(72)"7(a;2) — 

2 7 



(31) 



where 72 = 47me/(fcpmp) and 3:2 = fcpmpj/o/(47me). 
Finally the neutrino emissivity can be written as 



Qp-f-tpTr{Ei,, t) 



o 2 
-"p(73)"7(a;3)<Tp„-^ 



(32) 



Qp7^pee(7,*) = 2rpnp(7,t). 



(25) 



where 73 = 4Ev/{kpmpC^) and 3:3 = kpUipC^yo /{4Ev). 
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3.4- Synchrotron radiation 



3.5. Synchrotron self-absorption 



Classical synchrotron radiation is treated in several texts (e.g., 
Ginzburg &; Syrovatskii 1965). Using our normalisation, the 
rate rj^Xjf) at which a single electron electron emits photons 
into the frequency range da; is given by: 

,(,,^) = V!o^ i.[2./(36sin^7^)] (33) 

where b = B / Be denotes the magnetic field in units of the 
critical field Be = mlc^/{eh) = 4 • 414 x 10"G, at is the fine 
structure constant, 9 the particle's pitch angle and we have 
used the standard notation of Ginzburg and Syrovatski for the 
function F[x). Provided the energy of the emitted photon is 
much smaller than that of the emitting particle, the loss process 
can be considered continuous and represented by a first-order 
differential operator in the equation for electrons Eq. (11): 



d_ 
d-y 



"e(7,*) / 
Jo 



dx rj^Xjf) 



4 3 

-iB — [7'"e(7,<)] • 



(34) 



Here we have introduced the 'magnetic compactness': 

^==(^)-^ (^^) 

in which Ub = B^ /Stt is the magnetic energy density and have 
replaced sin^ 9 by its average (=2/3) for isotropic electrons. 
Equation (34) contributes an amount 4^B7min"e(7min, t)/3 to 
the source term Q^<'^°° for cool electrons. 

The source term in the photon equation can be found us- 
ing the '5-function' approximation, originally introduced by 
Hoyle (1960) in which the emission of a single electron is ap- 
proximated as monochromatic: 



■n{x,-y) w qoXo6{x - xo) , 



(36) 



and the quantities go and xo must be determined by imposing 
requirements on the accuracy of the approximation. To repro- 
duce the correct value for the total energy emitted by a single 
electron we find by integrating Eqs. (33) and (36) over x 



qoxo 



(37) 



One further constraint could be obtained by requiring the pro- 
duction rate of photons to be given precisely by the delta- 
function approximation. This leads to xo ~ 0-4667 • However, 
we prefer the simple expression: 



Xo = 67 



(38) 



which corresponds formally to requiring that the approxima- 
tion yields an accurate value for the moment J da;j;(a;,7)a;~'''^^. 
The appropriate photon source term follows simply from 
Eqs. (36) and (37): 



Ub 
36 
2 
3 



/OO 
d-yne{'y,t)6{x - 67^) 

iBb~^^^ x~^^^nc{'^ x/b, t) 



(39) 



Low energy electrons interact strongly with the radiation field 
if their Lorentz factor is such that the synchrotron photons 
they emit are reabsorbed within the source (McCray 1969, 
Ghisellini et al. 1988). The photon spectrum too is strongly af- 
fected below the self-absorption frequency (at which the optical 
depth equals unity). Such behaviour can be modelled using a 
second order derivative in momentum in the electron equation 
(Ghisellini et al. 1988). However, we do not wish to discuss in 
detail the heating and cooling processes of low energy electrons 
in this paper, nor are the details of the low energy spectrum 
important for our investigations. Consequently, we treat syn- 
chrotron self-absorption simply as a sink of soft photons. 

Starting from standard formulas for the absorption coef- 
ficient a„ (e-g-i Eq. (6.50) of Rybicki &; Lightman 1979), and 
using the 5-function approximation Eq. (36) one can readily 
derive 



Rciii,n~i[x, t) 

x~^^''b~^^''n.y[x, t) 

6af 



d_ 
d-y 



(40) 



7=(»,/6)l/2 



When only synchrotron emission and self-absorption are in- 
cluded, the photon kinetic equation (13) is: 



dn.y{^x^ i) 
dt 



+ 6 ^^'^x ^^''n^(x,t) 

6af 



d_ 
d-y 



(41) 



For an equilibrium electron distribution Jie oc exp(— 7/T) 
(where T is the temperature in units of meC^/fes) the sta- 
tionary photon distribution is the Rayleigh- Jeans distribution 
n^'' = 4iriZ(TT(Tnc//i)^Ta;, whereas for a power-law electron 
distribution one obtains a stationary self-absorbed spectrum 
n.y oc x^^'^ . These results agree with those found from an exact 
treatment of synchrotron emission. 

3.6. Inverse Compton scattering 

There is a close analogy between Compton scattering and syn- 
chrotron radiation, as has been pointed out by Felten &; Mor- 
rison (1966). In the case of synchrotron radiation, the classical 
treatment, in which the energy of the emitted photon is small 
compared to the energy of the electron, is adequate for our pur- 
poses. This is not so for Compton scattering, because events 
in the "Klein- Nishina regime" turn out to be important. We 
therefore divide our treatment of scattering into two parts. 

For collisions occurring in the classical or Thomson regime 
we use a method closely similar to that used for synchrotron 
radiation. The electron loss term is: 



A^s,t(7,*) = ^Ut-^ {y^ne{-y,t)) 



where 



f 

Jo 



(42) 



(43) 



is the energy density of photons with which the electron in- 
teracts via an inverse Compton scattering in the Thomson 
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regime i.e., photons of energy x' less than xi Ki 3/(47). Equa- 
tion (42) gives rise to an injection of cooled electrons at a rate 
4f7T7min"e(7min, *)/3- The photon source term in the Thomson 
regime is given by: 



^min[3/(4a;),3a;/4] 



Jo 



dx' 



x'~^'^x~^'^n^{^3x/4x',t)n^{x',t) (44) 



Photons with x > xt will interact with electrons of energy 
7 in Klein- Nishina scatterings. In order to treat this much more 
complicated case, we assume that an electron loses all its en- 
ergy in one such collision and joins the population of cooled 
particles. Approximating the cross-section by (t[x') ~ (Tt / x' 
where x' = -yx, the electron loss term can be written 



i-ics,KN(7i*) = / da; — 

J Xrr, 



(45) 



Because the scattered photons emerge with the energy of the 
incoming electron, the photon source term is simply 



Qics,KN(a;,*) 



p 00 

n^{x,t) I 



dx' [xx') ^n.y[x',t) 



(46) 



Note that these are only approximate expressions since the 
logarithmic energy dependence of the cross section has been 
neglected. Nonetheless, this does not affect our results because 
rates in the Klein-Nishina regime are in any case suppressed 
by the proportionality of the cross section to the inverse of the 
energy (the "Klein-Nishina" cut-off). 

3.7. Photon downs cattering on cold electrons 

Inverse Compton scattering describes the upscattering o{ pho- 
tons by relativistic electrons. In addition, cool electrons (7 Ri 
1), downs catterhaid photons. The Thomson optical depth of 
cool electrons is in our normalisation simply tt = N^°° . The 
rate of photon downscattering in the Thomson regime can then 
be found from the Kompaneets equation, assuming the tem- 
perature of cool electrons is zero and neglecting the stimulated 
scattering term: 



^cs(a;,*) = TT—[{Ax)n^{x,t)]. 
ox 



(47) 



(Ax) is the average energy shift of a photon with energy x 
colliding with an electron at rest. This quantity was calculated 
from the relation 



(Ax) 



/J;-"" dx'{x - a;')d<T(a;,a;')/dn 
/Jr" da;'d<T(a;,a;')/dn 



(48) 



where do'^x, a;')/dn is the differential cross section for Compton 
scattering (see, for example, Akhiezer &; Berestetskii 1969) and 
the limits of the integration x'^^^ = x and x'^^^ = x/[2x + 1) 
are determined from kinematical constraints. For x <C 1, 
(Ax) = x^/{2x + 1) ~ x^, which when substituted in Eq. (47) 
yields the usual Kompaneets equation (under the simplifica- 
tions assumed here). For a; < 1, (Aa;) departs from the a;^ 
dependence, as relativistic effects start becoming important. 
In order to take into account the spatial diffusion of photons 



inside the source, we follow Lightman &; Zdziarski (1987) and 
write for the photon escape time 



tje.c = 1 + -^f^^") ' 

where 

{1 a; < 0.1 

(l-a;)/0.9 0.1 < a; < 1 

a; > 1 . 



(49) 



(50) 



Therefore, when only the effects of photon downcomptoniza- 
tion are included, the photon kinetic equation (13) becomes 

dn~,(x,t) n~,(x,t) d ,,, . , , , 

+ -■ ^1 J( \ = ^X- Aa; a;,t . 51 

Ot 1 + -TTf(x) OX 

For a; > 1, Eq. (47) is no longer valid as the photons lose a 
large fraction of their energy in each collision. Compton scat- 
tering can then be considered as a catastrophic loss process and 
to deal with this case, we follow Svensson (1987) and write 

'Ccs(a;,*) = TTn^{x,t)/x. (52) 

The corresponding injection of (relativistic) electrons is 

QU-Y,t) = TTn^{T,t)h. (53) 

This injection comes from the upscattering of cooled electrons, 
and when 
Eq. (12). 



and when integrated over all energies contributes to z;^'*^""' in 



3.8. Photon-photon pair production 

This process acts as a sink of high energy photons, as well 
as an injection term of electrons. To calculate the loss rate of 
photons we follow Coppi and Blandford (1990): 
Photons of energy a; are lost at a rate 

/» 00 

Cl^~i^cc(x,t^ = n^(x,t^ I dx' n^(x' ,t^R^^(xx'^ (54) 
Jo 

where R^^ is a fit to the reaction rate given by 
- 1 

iZ^^(w) = • 652 —\n.{uj)H{ui - 1) (55) 

where H[y) is the Heaviside function. 

Photon-photon pair production is also responsible for the 
injection of relativistic electrons and positrons. Assuming, as 
in the case of Bethe-Heitler pair production, that the electron 
and positron emerge with equal energy 7, and noting that the 
photon-photon pair production process requires at least one 
hard photon of energy a; > 1, which interacts predominantly 
with a soft photon of energy around 1/a;, we find from conser- 
vation of energy 



27 = a; + 



(56) 



The injection term for electrons is then 

Q;^^ee(7,<) = 4n^(27,<) / dx' n.y{x',t)R.y.y{2-yx'). (57) 
Jo 
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3.9. Electron-positron pair annihilation 

This process is the inverse of photon-photon pair production 
described above. It acts as a sink of electrons and positrons and 
a source for photons. In this treatment we will consider only the 
"cosmic-ray" case discussed by Svensson (1982) and Aharonian 
et al. (1983) in which relativistic electrons/positrons annihilate 
only on the cool particles and not amongst themselves. With 
this assumption, the losses of electron-positron pairs can be 
written 

^ee^77(7l*) = TT^iZann (7)"e (Ti *) (58) 

where iZann is given by Coppi &; Blandford (1990) 

iZann = ^ [^-'I'+h.^] . (59) 

Photons produced by pair annihilation have {x) = (7 + l)/2 ~ 
7/2 so the photon term can be written 

QL^.yy(a;,t) = 27T:iZa„„(27)ne(27,t). (60) 

Finally the rate of Eq. (58) integrated over energy gives us 
the rate of cooled electron removal due to electron-positron an- 
nihilation, it contributes thus to £,'''^°° in Eq. (12). Since this 
rate does not include the annihilation of cooled pairs among 
themselves, we use an extra term 

r'"°'(t) = [Nt'"'\t)Y . (61) 

Furthermore, we assume that the rate of emission of the 
511 keV annihilation rate is given by the above equation, so 
we write 

QL^..(l,t) = liNr'it)]' . (62) 



be at 7min = 10°'^ . Electrons which cool through this boundary 
join the population of cool electrons. 

For the photons, the grid spacing in the logarithm of x 
is chosen to be twice as large as that in log 7 for protons (or 
electrons). Because the 5-function approximation is used in the 
synchrotron source term, Eq. (39), this grid spacing eliminates 
the need for interpolation; each electron grid point is associ- 
ated with a single photon grid point. Thus, since the softest 
photons are produced by electron synchrotron radiation, we 
take Smin = fcTmin- Howcvcr, inclusion of the pion producing 
processes requires the photon grid to be extended above the 
X corresponding to the synchrotron radiation of electrons of 
7max, because of the very hard photons produced in ir -decay. 
We therefore choose Smax = 7max and check that the photon 
bins adjacent to Smax remain practically empty throughout the 
run. 

Integrals over the photon and electron distributions are 
straightforwardly converted into summations using the trape- 
zoidal rule, but care must be taken in discretizing the first 
order derivatives with respect to 7. In the proton equation (6) 
and in the electron equation (11) these terms describe 'contin- 
uous' acceleration and losses and in order to avoid a numerical 
instability, an upstream difference must be taken. For particles 
undergoing losses 'upstream' means larger 7, whereas for the 
accelerating particles it means smaller 7. Thus, for electrons, 
which experience no acceleration, the value of the distribution 
at 7max is held constant at a negligibly small value. The value 
at 7min is computed. On the other hand, protons experience 
only acceleration at 7min, so that np(7min) is held constant 
and is related to the effective injection rate. At the upper end 
of the scale, losses always dominate for protons of 7max, so that, 
just as for electrons, np(7max) is held constant as a boundary 
condition. Similarly, in order to take into account the first or- 
der derivative with respect to the photon energy x (Eq. 47), we 
take an upstream difference while n^(a;max) is held constant at 
a negligibly small value. 



4. The numerical method 

4.1. Discretization 

When the physical process described in Sect. 3 are included, 
the system of kinetic equations (6), (11) and (13) is an 
integro-differential set for the distributions np(7,t), ne(7,t) 
and n^{x,t). Our approach to the numerical solution is to dis- 
cretize the variables 7 and x, and to integrate the resulting 
(stiff) set of coupled ordinary differential equations forwards 
in time. Because we use crude but computationally rapid ap- 
proximations to the physical processes we are able to integrate 
the equations using a standard NAG-library routine for stiff 
systems. 

The grid chosen for discretization is equally spaced in the 
logarithm of 7, the same grid being used for both protons and 
electrons. In order to accommodate a large dynamic range of 
the distributions np(7max,t), ie(7max,t) and n^{x,t) we use 
their logarithms in the numerical integration. In a typical run 
the resolution is 10 bins per decade for protons and electrons. 
This is a rather coarse grid; however, we found runs with higher 
resolutions to be excessively time consuming. The highest grid 
point 7max is adjusted so that the adjacent bins to ne(7max,t) 
and np(7max,t) remain negligibly small throughout the run. 
We find a value of 7max = 10 is adequate to ensure this for 
the models presented here. We choose the lowest grid point to 



4.2. Performance Checks 

In order to test both the code and the treatment of the physical 
processes, we have performed various runs for cases in which 
either an analytic solution is available, or in which we can com- 
pare our results with calculations in the literature. In the fol- 
lowing sections, when appropriate, we will refer to the electron 
injection compactness, defined in our normalisation by 

4 = ^^'^"""d7(7-l)Qe(7) (63) 

where 

Qe(7) = Qe,07~' for Tmin < 7 < Tmax /„ .x 

= otherwise 

is an external electron injection rate introduced for test pur- 
poses. 

4.2.1. Proton acceleration 

For the first test of the code we check the way the numerical 
method deals with proton acceleration. As stated above, pro- 
tons are effectively injected at a Lorentz factor 7inj = 10 ' . 
The solution of Eq. (1) in the case where = is given by 
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Eq. (2). Figure (1) shows the evolution of such a proton spec- 
trum in the case where tacc = tesc = 3tcross at times t = 10, 
20, 30 and 40 tcross (full lines). The dotted lines represent the 
analytical solution for this case. As can seen from this figure, 
the adopted numerical scheme is fairly accurate in tracing the 
time evolution of the cut-off' momentum. It also reproduces 
very precisely the expected power law, which has a slope —2. 
Some numerical diffusion is unavoidably introduced, which po- 




log(7) 



Fig. 1. The evolution of the proton spectrum in the no-loss case with 
tacc = tesc = 3 (full lines) at instants f = (a) 10, (b) 20, (c) 30 
and (d) 40 (all times are given in units of tcross). The dotted lines 
represent the corresponding analytical solutions. 

tentially inffuences the shape of any sharp features in the pho- 
ton spectrum. However, our acceleration model is in any case 
not accurate in the region of the spectral cut-off - a finite time- 
lag between scatterings as well as finite (but small) jumps in 
the particle energy (e.g., on crossing a shock) both lead to a 
smearing out of the sharp cut-off. 

4.2.2. Synchrotron radiation &; synchrotron self-absorption 

To check our treatment of synchrotron radiation we inject elec- 
trons with a power law spectrum at a constant rate given by 
Eq. (64) and include only the processes of synchrotron radia- 
tion, synchrotron self-absorption and photon escape. The re- 
sulting electron distribution can be found analytically (Karda- 
shev 1962). For s > 1 and 7 > 7min, there are three regimes. 
Defining the break point by 



7br(t) = 7max/(l +47max^Bt/3) 

we have: 

1. For 7 > 7max, 

ne(7,t) = for all t 



2. At early times i.e., for 7 < 7b 

3Qe,07' 



=(7,<) 



41b{s - 1) 
which, for UBlt/Z < 1 



7 7 



the spectrum is 
4^B7*/3)''"^ 



(67) 



is a power law of index s in 7 
with an amplitude which increases linearly with time: Ki 

3. For 7 > 7br, electrons have had time to cool, and the spec- 
trum steepens: 

3Qe,0 



=(7,<) 



-{s-l) 



-{'- 

Tmax 



(68) 
cal- 



41b{s - 1)72 ^' 

Figures (2A) and (2B) show the numerically 
culated electron and photon spectra at times t = 
10~ , 10~ , 10~ , 10~ and 10 (in units of tcross), choosing 
b = 10"% is = 0.075, s = 2, 7inin = 3 and 7inax = 3.10\ 
The electron spectra show breaks at energies which agree well 
with Eq. (65). For 7 <C 7br the electron distribution is approx- 
imately a power-law Jie oc 7~", has the same index as the in- 
jection term (a = s) and increases linearly with time, whereas 
for 7 ^ 7br electrons have cooled giving a time-independent 
spectrum of index a = s + 1 = 3. 

The photon spectrum shows similar behaviour. Here, how- 
ever, there are two breaks. The one at high energies which 
we shall call Sbr is related to 7br by Sbr = 67bri that at low 
energies, called Sssa, is caused by synchrotron self-absorption. 
Photons with Sssa <C a; <C a^br have a power-law distribution 
n-y oc with p = 1.5, in agreement with the standard for- 
mula p = (a + l)/2. Photons with x > Xhi correspond to the 
'cooled' part of the electron distribution (7 > 7br) and have 
a power-law of index of p = 2 again in agreement with the 
standard formula (in this case with a = s + 1 = 3). Photons 
with X < Sssa show the characteristic self-absorbed distribution 
n-y oc x^^'^ which, as can be seen from Eq. (41), is produced by 
any power-law electron distribution, independent of the value 
of s. It is interesting to note that while Xhi decreases with time 
(since 7br also decreases), Sssa increases. This happens because 
as the system evolves in time, the number density of low en- 
ergy electrons and, hence, the optical depth to self-absorption 
increases. 

4.2.3. Inverse Compton Scattering 

A similar approach to the above can be taken for the process 
of inverse Compton scattering. We take a constant electron 
injection given by Eq. (64) and assume the scattering is with 
photons of an external black-body field of temperature = 
fcT/meC and compactness ^bb. The external photon source 
can be written 



"bb 



454b 



^404 g»,/G _ J- 



(69) 



We then solve numerically the kinetic equation for electrons 
retaining only the Compton loss terms, and the kinetic equa- 
tion for photons including inverse Compton scattering and es- 
cape from the boundary. In order to illustrate the effects of 
the transition from the Thomson to the Klein- Nishina regimes 

(65) 

(which, in our case, are described by two different electron loss 
and photon gain rates) we solve the equations keeping all pa- 
rameters the same and increasing only the upper limit of the 
electron injection 7max. In the particular case considered here 
we choose = 10~ and 4b = 1/3 for the parameters of the 
(^^) external photon field, and s = 2, 7min = 3 and Qe,o = 10~^ for 
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Fig. 2. (A) the electron distribution for synchrotron cooUng and con- 
stant power-law injection. (B) the corresponding photon distribution 
including escape and synchrotron self-absorption. Here b = 10~^, 
= 0.075, 5 = 2, 7min = 3 and 7max = 3.10^. The curves show the 
evolution of the electron and photon distributions towards steady 
state. Curve (a) is t = 10"'* after electron injection, (b) is for 
t = 10-3, ig fQj. ^ ^ 10-^, (d) is for t = 10-1 ig fQj. 

t = 10. All times are expressed in units of tcross* 




log(7) 
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Fig. 3. Electron steady state and photon spectra in the case where the 
electron cooling is provided by inverse Compton scattering on ther- 
mal photons of temperature = 10—^ and compactness = 1/3. 
The electrons were assumed to be injected with a power law of slope 
s=2 between 7min = 3 and 7max= (a) 3.10^, (b) 3.10^, (c) 3.10^, 
(d) 3.10^ and (e) 3.10^. 
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the electron injection parameters. The results are presented in 
Fig. (3)for7inax of (a) 3.10^ (b) 3.10\ (c) 3.10^ (d) 3.10^ and 
(e) 3.10^. Since the blackbody photons can be approximated 
as having an energy (sbb) — 30 = 3.10"^, we have that for 
cases (a) and (b) all scatterings occur in the Thomson regime, 
whereas for the other cases, the more energetic electrons scat- 
ter in the Klein Nishina regime. Figure (3A) shows the steady 
state electron distribution and Fig. (3B) the steady state pho- 
ton distribution. For cases (a) and (b) one obtains, just as in 
the case of synchrotron cooling, a power-law electron distrib- 
ution of index a = s + 1 = 3 and, correspondingly, a photon 
spectrum of index p = 2. For the other cases scattering in the 
Klein-Nishina regime is important and the electron spectrum 
becomes harder for 7 > 10^. In accordance with previous work 
(Jones 1994, Blumenthal and Gould 1970, Blumenthal 1971, 
Zdziarski 1989) we find an index a~s — l = las 7max becomes 
large. However, as can be seen from Fig. (3B), the photon dis- 
tribution does not break but retains itsp = 2 slope throughout. 
This is to be expected for the special case of s = 2, because 
in the Klein Nishina case the photon distribution has a slope 
p = a + 1 = 2, i.e., the radiated photons have essentially 
the injected electron index (Blumenthal &; Gould 1970) and in 
the Thomson regime one has p = (a + l)/2 = 2. Although 
our approximate treatment reproduces the above properties 
quite well, the abrupt switch we use between the two scatter- 
ing regimes produces artificially sharp breaks, such as can be 
seen in Fig. (3A). 



4.2.4. Electromagnetic cascades 

As a test for Compton scattering, photon-photon pair produc- 
tion, electron-positron annihilation and Compton downscatter- 
ing, we present a series of runs showing pair cascades, which 
have been investigated by many authors [see Svensson (1989) 
for a review]. Electrons injected with s = 1.5, 7min = 3 and 
7max = 7.10^ are allowed to cool on an external photon field 
with = 10~ and ^bb = Is- For electron compactnesses 1^ of 
(a) l/(4ir), (b) 30/(4ir) and (c) 1000/(4ir) we ran the code until 
a steady state was reached. The electrons cool on the external 
photons as in the case examined in Sect. (4.2.3). However, es- 
pecially when the compactness is high, the resulting 7-rays 
do not escape freely. Instead, they produce electron-positron 
pairs upon interacting with soft photons. These pairs, in their 
turn, produce more 7-rays, initiating an electromagnetic cas- 
cade. The computed photon spectra are shown in Fig. (4). In 
case (a), where the compactness is quite low, the computed 
spectrum is the same as that expected in the pure inverse 
Compton scattering case, i.e., there is no pair feedback. For 
intermediate compactness (case b) there is some pair produc- 
tion and the spectrum steepens because of the redistribution 
of the luminosity to lower energies by pair production. Finally 
for high compactness (case c) there is substantial pair produc- 
tion and 7-rays of a; > 1 are effectively absorbed. However, the 
steepening of the spectrum at energies a; < 1 is due not to pair 
production, but to the downscattering of photons on cooled 
pairs - an effect which is important for l/rrf < x < 1. We 
have compared our results with those found by Lightman and 
Zdziarski (1987), for the same initial conditions (their Fig. 3 - 
note that there is a 4ir difference in the definition of compact- 
ness between Lightman and Zdziarski and the present work). 
Overall we found a very good agreement, the only major dif- 
ference being in case (c) around x ~ 10~^. This effect can 



o 



1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 




(c) 




(b) \ 


/, /l /, , , 


fa) ^ W 

, , , 1 , , , 1 , , , 1 , , w 



log(x) 



Fig. 4. Photon spectra obtained from electromagnetic cascades. The 
injected electrons have a slope s = 1.5 between 7min = 3 and 
7max = 7.10^ and they cool on an external photon field with 
= 10~^ and ^bb = ^e. Curve (a) corresponds to an electron com- 
pactness of l/(47r), curve (b) to = 30/(47r) and curve (c) to 
4 = 1000/(47r). 



be attributed to the upscattering of soft photons by thermal 
electrons - a process we neglect. 

5. The Pair-Production— Synchrotron Instability 

We turn now to an example of the use of the full code. The 
system of the kinetic equations (1), (11) &; (13), is subject 
to various feedback effects and in this section we present one 
of them, namely the 'pair-production synchrotron instability' 
(PPS). This feedback effect manifests itself when the number 
density of relativistic protons with energy above a threshold 
[7 > 7crit = (2/6)^''^] becomes larger than a certain critical 
value. In this case the relativistic protons produce pairs which 
subsequently radiate synchrotron photons energetic enough to 
provide the next generation of targets for the relativistic pro- 
tons. This eventually leads to a runaway - (KM92). 
The input parameters for the simulation are 

1. the strength of the magnetic field, expressed as the ra- 
tio b = B/Bc. Although the threshold 7crit depends only 
weakly on 6, the growth rate of the instability is more 
strongly affected. In the simulation, we choose for conve- 
nience a value of 4.4 x 10^ G (6 = 10~^), rather higher than 
those normally estimated for AGNs. 

2. the magnetic compactness is, which, given the mag- 
netic field, determines the size of the emission region 
via Eq. (35). For this particular run we choose a size of 
1.5 X 10^^ cm, corresponding to a light crossing time of 
about one hour. Such a value is appropriate for a rapidly 
variable Seyfert galaxy. 
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3. 



the acceleration time tacc- The pair-production synchrotron 
instability threshold can be derived assuming acceleration 
is a slower process than all others (except escape, which 
has exactly the same timescale in our case, leading to a 
jip oc 7~ spectrum). Our treatment becomes inconsistent 
when the acceleration time is less than a light crossing time, 
since this implies small regions within the source where 
acceleration and losses compete, whereas we treat each of 
these processes as spatially homogeneous. Thus, we choose 
an acceleration time tacc equal to three light crossing times, 
the compactness of the total nonthermal luminosity ^tot, 
which is related to the quantity no defined in Eq. (8) by 
the relation (9). This parameter determines whether or not 
instability occurs. From KM92, we have for the threshold 
value Jicrit : 

-1 



1 



,1/3 



f 

Jo 



-5/3 



(70) 



where Upe is the pair production cross section (in units 
of (Tt). For no > ncrit there exists a critical Lorentz fac- 
tor 7max for the upper cut-off' of the proton distribution, 
above which the system goes unstable. In the simulation 
we choose no = l.lncrit and the instability takes off soon 
after the the maximum Lorentz factor exceeds the value 
Tcrit = = 1.3 X 10\ 

5. the background density of protons. We assume here that 
proton acceleration takes place in an environment where 
the thermal plasma density can be considered negligible, 
as may be the case in an accretion disk corona. Thus the 
only proton-proton interactions we consider are those of 
the relativistic protons amongst themselves. 

Fig. (5) shows the evolution of the photon compactness 
(solid line) and also a quantity dp which we call the 'proton 
column-depth' and which gives a dimensionless measure of the 
energy contained in relativistic protons in the source: 



dp 



3m, 



f 



d7 77ip 



(71) 



Initially, a few seed photons are produced by proton-proton 
collisions (i.e., in 7-rays from neutral pions and synchrotron 
photons from the electrons/positrons injected from the decay 
of charged pions - cf. Section 3.1). However, at about 20 cross- 
ing times the upper cut-off of the proton distribution becomes 
larger than 7crit. The PPS instability then takes over and re- 
sults in a rapid increase in the photon density at about 25 
crossing times. The protons react to this growth after only a 
few crossing times, their energy density being reduced as a 
result of cooling on the produced photons. Because the pairs 
produced in this process cool rapidly, the photon density de- 
creases on a time scale of about the light crossing time. This 
enables the acceleration process to act again, increasing the 
proton 'column-depth' until pair production restores the pho- 
ton density. The resulting oscillations are damped for these 
parameters, and the the situation reaches a stationary state 
after several crossing times. The photon compactness in the 
stationary state is only roughly 50% of the total compactness 
^tot (shown by the dashed line). The remaining nonthermal 
power escapes the system in the form of either neutrinos or 
neutrons. The interrelation between i-y and dp is displayed in 
Fig. (6). 

Figure (TA) shows the proton spectrum at five values of t 
chosen in such a way as to show the onset and saturation of 




Fig. 5. Photon compactness (full line) and proton column-depth (dot- 
ted line) versus time expressed in units of crossing time in the case 
where the pair/synchrotron instability operates. Injection parame- 
ters are no = l.lncrit) *acc = Stcross and b = 10~^. The total com- 
pactness ^tot [see Eq. (10)] is shown by the dashed line. 




Fig. 6. Photon compactness versus proton column-depth for the in- 
jection parameters specified in Figure 5 and in the text. 
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the instability. These snapshots are for ti = 23, 28, 33, 38 and 
43 (i = 1 ... 5) crossing times after the run was started. The 
corresponding photon distribution is shown in Fig. (7B). The 
photon compactness at these five instants is i-y ~2, 4, 13, 28 
and 21 respectively. At ti (at which point the photons have 
still a low compactness - i-y ~ 2, and have not yet become 
important for proton losses) the protons have a power law dis- 
tribution up to 7max (~ 2 X 10 in this particular case), i.e. at 
this stage their upper cut-off' has just moved at values higher 
than 7crit(— 10 ). The proton slope is —2.02, a little steeper 
than the canonical value of —2 and this can be attributed to 
the effects of proton-proton losses - cf. Section 3.1. At t2 and 
ti the proton energy continues to increase with the photon 
energy density also increasing but still too low for effective 
proton losses. However, as the photon energy density increases 
(for t > ti) losses become important and the protons cool down 
to values close to 3.10 , not too far from their final steady state 
distribution. 

On the other hand the photon spectrum has the following 
characteristics: As the instability manifests itself the photons 
have a spectral index of —1.8, close to the value of —5/3 pre- 
dicted in the linear theory - see KM92. At high frequencies it 
shows a turnover which is given approximately by the relation 
Sbr — fcTmax ~ See scctious (3.2) and (3.4). As 7max moves to 
higher energies (i.e., for ti to ts), so does Xhi- However, as the 
protons cool (for t = ti), Xhi moves to lower energies. We find 
that this break varies between 50 - 200 keV as it is to be ex- 
pected from the formula for Xhi used above. At low frequencies, 
the photon spectrum shows strong synchrotron self-absorption 
as expected from such a compact source. Furthermore, as the 
source flares up, the self-absorption frequency increases. Again 
this is to be expected since during the linear stage of the insta- 
bility the number density of electrons in the system increases. 

Fig (8) depicts the photon number spectral index a over the 
range 2-10 keV as a function of the photon compactness i-y. The 
slope is a little steeper than the analytically calculated slope of 
— 5/3 (KM92). It is worth noting that the spectrum becomes 
steeper as the photon luminosity decreases after reaching the 
peak. 

The dimensionless neutrino emissivity Q„ once steady state 
has been established is shown in Fig. (9). This is calculated us- 
ing Eqs. (19) and (32). The neutrinos produced in pp collisions 
dominate at relatively lower energies (£^„ < 100 GeV) and have 
a slope of ~ —2, reflecting the parent proton distribution. The 
neutrinos produced in p7 collisions dominate at higher ener- 
gies, however they fall abruptly above a few TeV. 

In this example the evolution of a system follows closely 
that predicted for the PPS instability. Out of all the processes 
described in Section 3 the leading ones are proton-photon 
pair production and synchrotron radiation. The effects of 
other processes such as photopion production, inverse Comp- 
ton scattering and photon-photon absorption remain negligible 
throughout. 

6. Summary 

In this paper we develop a technique suitable for the inves- 
tigation of AGN models in which the basic source of lumi- 
nosity is the acceleration of nonthermal protons. Once these 
have reached sufficiently high energy, they are responsible for 
the injection of electron/positron pairs and photons. Such 
models have been suggested some time ago (Protheroe and 
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Fig. 7. (A) The proton distribution Tip at five different times for 
the injection parameters specified in Figure 5. The soUd Une is for 
t = 23, the dotted Une for t = 28, the short-dashed Une for t = 33, 
the long-dashed Une for t = 38 and the dot-dashed Une for t = 43. 
AU times are expressed in units of tcross* (B) The corresponding 
photon spectrum. 
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Fig. 8. The 2-10 keV photon number index a versus photon compact- 
ness for the parameters given in Figure 5 and in the text. 
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Fig. 9. The dimensionless neutrino emissivity Qv as a function of 
the neutrino energy in GeV once steady state has been established. 
The dotted line corresponds to the emissivity due to proton-proton 
collisions, the dashed line is the emissivity due to proton-photon 
collisions and the full line is the total emissivity. The injection para- 
meters are as given in Figure 5 and in the text. The true emissivity 
(number of neutrinos emitted per second per unit energy per unit 
source volume) is ai, = Q v / {rn^cB? a^^^ 



Kazanas 1983, Kazanas and Ellison 1986), but the new aspect 
of our work is the self-consistent treatment of the interactions 
between the three important components of the model: pro- 
tons, electrons/positrons and photons. 

To achieve this, we follow the kinetic equation approach us- 
ing spatially averaged distributions and an 'escape probability' 
formalism. Self-consistency is achieved by including source and 
loss terms in these equations which account for all of the im- 
portant interactions between the components. Furthermore, we 
avoid prescribing the injection of nonthermal power by using an 
explicit model for the acceleration of the protons. This results 
in a rather constrained parameter space for our problem. The 
system is fully specified when the details of proton acceleration 
and the injection of low energy seed particles are determined. 
The electron and photon equations contain no additional free 
parameters. Thus, the external parameters usually employed 
in pair plasma problems, such as the energy with which non- 
thermal pairs or photons are injected and their compactnesses 
(see Lightman and Zdziarski 1987 and Svensson 1987) are in 
our case absent. 

The resulting system is sufficiently complicated that a num- 
ber of checks of both the physically motivated approximations 
used for the elementary processes, as well as the numerical 
method, is essential. These checks indicate that the method 
provides a useful tool for more detailed applications to AGN. 
As well as discussing tests, we present in this paper one exam- 
ple of results which are obtained with all processes included 
- one in which the feedback loop of the the pair /synchrotron 
instability is expected to operate. We find this instability is in- 
deed present, and is not significantly suppressed by processes 
such as proton-proton interaction or photo-pion production. 
The prediction of an X-ray spectral index of —5/3 which stems 
from a linear analysis of this instability (KM92) is shown to 
be unaffected by the additional processes we include, and to 
extend well into the nonlinear phase. 

Another interesting feature of the pair /synchrotron insta- 
bility is that the photon spectrum breaks at energies between 
50 and 500 keV and this agrees generally well with the re- 
cent OSSE observations of Seyferts (see Maisack et al. 1993). 
While there are certainly other models to explain this feature 
(see Zdziarski et al. 1993, Titarchuk &; Mastichiadis 1994), the 
fact that, on the one hand, we can explain such a feature in 
a simple way within the framework of our model and, on the 
other, reproduce the correct 2-10 keV spectral slope, makes the 
present approach promising. 

A limitation of the method is that it is not capable of treat- 
ing the details of radiation transfer through the source, using 
instead an escape probability formalism. This problem is com- 
mon also to other models treating pair plasmas in the kinetic 
equation approach. Furthermore, the microphysical processes 
were treated using simple approximations and some processes 
such as Coulomb collisions and bremsstrahlung were neglected 
altogether. However, these processes are important only in a 
rather limited region of parameter space (see Svensson 1982 
and Coppi 1992), which we avoid. 

Other effects we do not try to include in this work include 
the possible radiation from protons escaping downstream and 
the effects such radiation might have on the acceleration zone. 
This is justified if, for example, escaping particles are accreted 
into the black hole. Finally, we do not attempt to treat the 
effects of neutrons produced in hadronic collisions and their 
possible feedback on the system. However, at least for the para- 
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meter space where the pair /synchrotron instability dominates, 
protons lose the bulk of their energy in proton-photon pair 
production and neutrons are, therefore, unimportant for the 
evolution of the system. 
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